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Abstract 

The Luria-Delbriick mutation model has a long history and has been mathemat- 
ically formulated in several different ways. Here we tackle the problem in the case 
of a continuous distribution using some mathematical tools from nonlinear statistical 
physics. Starting from the classical formulations we derive the corresponding differ- 
ential models and show that under a suitable mean field scaling they correspond to 
generalized Fokker-Planck equations for the mutants distribution whose solutions are 
given by the corresponding Luria-Delbriick distribution. Numerical results confirming 
the theoretical analysis are also presented. 

Keywords: Luria-Delbriick distribution, kinetic models, mutation rates, Fokker- 
Planck equations 

1 Introduction 

Application of tools from nonlinear statistical physics to describe different biological phe- 
nomena started to gain popularity in the recent years and it represents one of the major 
challenges in contemporary mathematical modeling [H [5| [TT| I22j. In particular, powerful 
methods borrowed from kinetic theory of rarefied gases have been successfully employed 
to construct master equations of Boltzmann type, usually referred to as kinetic equations, 
describing the emergency of universal behaviors through their equilibria [71 [23l [26] . 

An important example of emergent behavior describes building of tumors by cancer 
cells and their migration through the tissues \IU[ fTTj \T2\ . Another famous example to 
consider in this contest is the classical Luria-Delbriick mutation problem |19l [3T] . Cer- 
tainly, the basic entities in these examples differ from the physical particles in that they 
already have an intermediate complexity themselves. However, for some specific phenom- 
ena, the statistical behavior of the system is related to the peculiar way entities interact 
and not to their internal complex structure. 

First experimental and theoretical analysis of mutation process in bacteria was pub- 
lished by Luria and Delbriick in 1943 [19] (Nobel Prize in Physiology and Medicine, 1969). 



* Applied Mathematics Department, Tel Aviv University, Israel (ekashdaii@post.tau.ac.il) 
^Mathematics Department, University of Ferrara, Italy (lorenzo.pareschi@unife.it). 



1 



The goal of the study conducted by Luria and Delbriick was to estimate the mutation 
rate in bacterial populations by observing the fraction of the final cells that carry a muta- 
tion. In their experiment, often called fluctuation test, they have shown that in bacteria, 
genetic mutations arise in the absence of selection, rather than being a response to selec- 
tion. They explained these results mathematically by the occurrence of a constant rate of 
random mutations in each generation of bacteria growing in the initial culture. 

The distribution of the number of mutants in a culture that has been grown under 
conditions in which these mutations did not confer a selective advantage to the cells 
bearing them, came to be known as the Luria- Delbriick distribution. 

In multicellular organisms, connection between mutagenesis and carcinogenesis is broadly 
accepted (see, for instance, [24]) and the Luria-Delbriick distribution plays an important 
role in the study of cancer, because tumor progression depends on how heritable changes 
(mutations) accumulate in cell lineages. In particular, the Luria-Delbriick model focuses 
on the distribution of mutations in an exponentially expanding clone of cells. We should 
mention here that several models describing the Luria-Delbriick experiment have been pro- 
posed. The most famous is the Lea-Coulson formulation ^18j based on the introduction of a 
random rate of mutations instead of constant rate implemented in [19]. Other formulations 
have been introduced in [21 133] . We refer the reader to [31] for a comprehensive review 
of the subject. The mathematics behind the Luria-Delbriick experiment (most often in 
the Lea-Coulson formulation) has been studied by several authors [H [T2 | [T5 | \T6 [ \20 [ [32]. 

The specific shape of Luria-Delbriick distribution is strictly connected to the partic- 
ular formulation of the mutation dynamics. In general, no analytic expressions of such 
distributions are available, however, exact expressions may be available for the character- 
istic functions, the probability generating functions and the moments. This indeed is the 
case of both the original (Luria-Delbriick) and the Lea-Coulson formulations considered 
in this manuscript. For such a reason a lot of effort has been invested in the development 
of computational methods and asymptotic approximations of Luria-Delbriick distribution 
(see ^31j and the references therein). 

Here we introduce kinetic master equations, in the spirit of [26j, where the number of 
mutants is considered as a continuous variable instead of a discrete one. The behavior of 
the resulting differential models are then studied by using a suitable mean field asymptotic 
scaling, which permits to recover the exact moments of the original formulations. As the 
limit we obtain generalized Fokker-Planck equations and, for the classical Luria-Delbriick 
setting, show that the long time behavior is characterized by the corresponding Luria- 
Delbriick distribution. The same evidence is give numerically also for the Lea-Coulson 
setting. 

The rest of the paper is organized as follows. In the next section we consider the stan- 
dard Luria-Delbriick experimental settings and present the corresponding kinetic master 
equations. Then we introduce a mean field scaling that permits to derive a generalized 
Fokker-Planck model whose solution is given by the original Luria-Delbriick distribution 
of mutants. In section 3 we consider the case of a variable mutation rate as proposed by 
the Lea-Coulson formulation. Derivation of the generalized Fokker-Planck model for this 
case is also presented and its solution is discussed. Simple diff^usion approximations of 
the mutation process are also derived in this section. Then, in section 4, we present nu- 
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merical experiments confirming our theoretical analysis. Concluding remarks and possible 
extensions of this work are reported in the last section. 



2 The Luria-Delbriick formulation 

The mathematical formulation of the Luria-Delbriick model is based on the following 
assumptions 

• The process starts at time t = with a number A'^o of normal cells and no mutants. 

• Normal cells replicate at a constant rate a. 

• Mutation occur randomly at a rate characterized by a Poisson process with intensity 
proportional to the total number of normal cells. We denote with /j, the mutation 
rate per cell. 

• The offspring of every mutant cell replicate at a constant rate /3. 

As a consequence of the above assumptions if we denote by N(t) the number of normal 
cells we obtain the differential problem 

^^ = aN{t)-fiN{t), N{0) = No, (1) 

which gives 

N{t) = A^oe^""'')*, t > 0. (2) 

The fundamental question in the above model is the estimation of the distribution of 
the number of mutants in time. Such distribution, in fact, is of vital importance when 
estimating the mutation rates. We shall denote by f{m, t) the probability density of having 
m mutant cells at time t. 



2.1 A kinetic model for mutations 

In the sequel we will assume that the random variable m characterizing the number of 
mutants at time t takes values in R+. Such an assumption, as we will see later, permits us 
to tackle the problem with some tools borrowed from kinetic theory [6l [71 [23| 126) . Thus 
we have 

/ /(m, t) dm = 1, 
and denote by M{t) the expected number of mutant cells at time t 

M{t)= I f{m,t)mdm. (3) 
Jr+ 

We consider the following microscopic dynamic for the random variable m 

m' = m + /3m + rj, (4) 
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where ?? > (absence of backward mutations) is a discrete random variable distributed 
accordingly to a Poisson density p{N{t),rj) with mean fj,N(t). 

By standard arguments we can write the following kinetic master equation for the time 
evolution of the distribution of mutant cells 



df{m,t) 
dt 



E 

i=0 



B,m^^{N{t),7]i)jf{'m, t) - B^^m'{N{t),rji)f{m, t) 



(5) 



where 'm = {m — r]i)/ J , J = 1 + /3, and the transition rates are given by 

B,^^^{N{t),r^i) = p{N{t),r]i)xi!m > 0), Bm^m'{N{t),r]^) = p{N{t),7],i)x{m' > 0), 

with xi^) the indicator function of the set /. The above equation is complemented with 
the initial data 

/(m,0) = /o(m). (6) 

Since the random variable r/j takes values only in M+ the weak form of the kinetic 
equation yields the simpler representation 



d_ 
di 



f{m,t)ip{m) dm = ''^^p{N{t),7]i) / f{m,t)[ip{m') — ip{m)]d'm 



i=0 



(7) 



where ip(m) is a smooth function. 

It is easy to see that taking ip(m) = 1 we get the consistency condition 

- fim,t)dm = 0. (8) 

The choice ip(m) = m yields the evolution of the mean number of mutant cells and gives 

f{m,t)mdm = ''^^p{N{t),7]i) / f{m,t)[l3m+r]i]dm = (3 / f{m,t)mdm—fj,N{t), 



d_ 
dt 



i=0 



which corresponds to the simple ODE 



^^ = /3M(t) + /xiVoe("-'^)*. 
dt 



(9) 



The above equation can be solved exactly and, starting with initial data M(0) = Mq, gives 

/3 = a — fi. 



M{t) 



a — fj. — /3 



(10) 



Note that the above expression for the mean coincides with what found in the literature 
for the Luria-Delbriick distribution 
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Similarly we can compute the evolution of the variance. Starting from ip{m) = m? we 
have for the second moment M2{t) the equation 



9_ 
di 



p CO „ 

/ f{m,t)m^dm = Vp(iV(t), r/i) / f{m,t)[^ {2 + ^)m'^ + ij^ + 2{1 + ^)mr]i] dm, 
Jr+ ~^ Jr+ 

= /3(2 + /3) / f{m,t)m^dm 
Jr+ 

+ <r]^ > +2{1 + (3) fiN{t) f{m,t)mdm, 

Jr+ 



where < r/^ >= fiN{t){l + fiN{t)) denotes the second moment of r/. 
Thus the variance 

V{t)= f f{m,t)m^dm- M{tf, 
Jr+ 

follows the differential equation 
dV{t) 



dt 



2^V{t) + rM2(t) + fiN{t){l + 2PM{t) + fiN{t)), (11) 



which differs from the variance of the Luria-Delbriick distribution. The exact analytical 
solution to (jlip can be easily computed, we omit it for brevity. 

To analyze the long time behavior let us consider the case Nq = 1, Mq = 0, which 
represents the situation, where the mean numbers of normal cells and mutant cells are 
initially one and zero respectively (a standard Luria-Delbriick setup [19j). Note that, 
when the growth rates are identical a = /?, the average total number of cells is B{t) = 
M{t) + N{t) = e^* which is simply a birth process at rate {3. 

Let us define next the concentration of the mutant cells as p{t) = M{t) / B{t). The 
long time behavior of p according to ()10p is given by 

!1, I3>a-p, 
o (12) 
a-/3 

The behavior is sketched in Figure 1. Mutants become predominant when f3 > a — p. 
It is interesting to revise that the asymptotic behavior of p shows exponential convergence 
for all values of /? except for (3 = a — p, where we have 



l + pt 



For 13 < a — p mutant cells represent a fraction p/{a — (3) of the population and no 
predominant behavior is observed. 

2.2 Derivation of the generalized Fokker-Planck model 

In order to study the asymptotic behavior of the model we introduce the following scaling. 
We set 

f(m, t) = f{m, t), T = et. (13) 
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Figure 1: Long time behavior of mutant cell concentration as function of time 



We are interested in the behavior for small values of e such that 

lim — = 7, lim — = 71 lim — = v, (14) 

e,/3^0 e £,/3i-)-0 e e,At^O £ 

where /3i = a — 

Roughly speaking the above scaling limit implies m' m and corresponds to the realis- 
tic assumption of considering the the birth-mutation dynamic as a slow process originated 
by many small variations of the number of mutants. For this reason we will refer to the 
limiting dynamic as a mean field dynamic. In particular, as we will show, it is consistent 
with the evolution of the variance of the original Luria-Delbriick model. 

To this aim, let us first point out that such scaling does not affect the evolution of the 
time rescaled mean M(r) which in the limit ()14p is governed by 

^MSZI = ^M{r) + z^iVoe^i". (15) 
ar 

For the time rescaled variance V{t) in the limit e — )• we obtain the simpler equation 

^Kill = 2-iV{T) + vN^e^^^ (16) 
ar 

Clearly the above equation can be solved exactly, and for F(0) = gives 



V{t) 



^ViVoe'^^(e(^^-2^)^ - 1), 27/71, 

71 - 27 (17) 

IiNqtc'^^, 27 = 71, 
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which is the same expression of the variance for the Luria-Delbriick distribution |31j . 
The time scaled distribution f{m, r) satisfies the equation 



— / f{m,T)(p{m)dm = -y^p{N{T),rii) [ f{m,T)[ip{m') - (p{m)]dm. (18) 
Let us expand (p{m') in Taylor series 



fc=i 



and substituting into ([TH]) we have 



o 



f{m,T)(p{m) 



" k=l "' i=0 



dr 

Using the binomial formula we can write 



d^ = -y2ny2pi^i^)^Vi) / f{m,T){l3m + 7]i)''^^''\m) dm. (19) 



{/3m + rji)'' = E (^) /^'"i'^f"' 

n V J / 



i=o 

and then 

00 ^ A; 



— — /" f (m, r)(/3(m) dm = — I | — — [ f(m,T)m^(p'^^\m)dm, 



k=l j=0 

where < r)^~^ > denotes the (A; — j)-th moment of rj. 

We are interested in taking the limit e, /3i, /3, /i — > 0. Now since r]i is given by a Poisson 
process 

00 

and thus 

< v^^^ > 

lim — '- = i^N{t), yj <k. 

e,p.^O e 

In the limit we finally obtain 



— / f{m,T)ip{m) dm — ^ [ mf{m,T)ip'{m)dm = 'S^—N{T) [ f{m,T)Lp^^\m) dm, 
OT Jk+ ^ fe! Jbl+ 

which corresponds to the generalized Fokker-Planck equation 

—f{m, t) + 7— (m/(m, r)) = z.iV(r)/:i,M(/(m, r)), (20) 



where 



^ ( ■\\k P)(k) 

CKM{hm,r)) ■.= Y.^-^d^)hm,r) (21) 

k=l 

is the Kramers-Moyal operator. 

Let us scale our solution accordingly with 

g{m, t) = u{T)f{mu{T),T), u{t) = (22) 

Then g{m, r) satisfies the equation 

-g{m, r) = uN{r) J] ^^^^ d^A"^^ (^3) 

fc=i 

Now, we introduce the Fourier transform 

m.r)= I gim,T)e-'"'^dm, 

JR 

to obtain 



o oo / -] \ k 



■yT\k 

k=l 

i/iVoe^5(C,r)(e-^«^"^^-l) 



The solution of the above differential equation can be written in the form 

m r) = m) exp (^uNo £ (e'^^^""^ - l) dz^ , (24) 

with 

50 (0= [ go{m)e-'^^ dm. 

Note that equation (|24|) is exactly the characteristic function of the Luria-Delbriick distri- 
bution [31]. This shows that the solution of the kinetic model d?]), in the asymptotic limit 
described by equations (fTB|) -(fT ^ . coincides with the classical Luria-Delbriick distribution 
characterized by (p^ . 



3 The Lea-Coulson setting 

The original formulation by Luria and Delbriick assumed a deterministic growth rate for 
mutant cells, which seemed too stringent an assumption to allow for efficient extraction 
of reliable information about mutation rates from experimental data j^Slj. A slightly 
different mathematical formulation was proposed by Lea and Coulson [18], who adopted a 
preferential attachment process (or Yule process) with birth rate /3 to describe the growth 
of mutant cells. Today the Lea-Coulson formulation is the prominent model for the study 
of mutation rates and is now commonly referred to as the Luria-Delbriick model. 
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3.1 The kinetic equation 

The formulation of the Lea-Coulson assumptions is based on considering the microscopic 
dynamic 

m' = m + ■!? + ?7, (25) 

where now i9 > is a discrete random variable distributed accordingly to a Poisson density 
p{m, t?) with intensity /3m. 

The weak-form of the kinetic equation now reads 

— / f{m,t)ip{m)dm = p{N{t),r]i) p{m,'&j)f{m,t)[ip{m')-(f{m)]dm, (26) 
ot M+ Jr+ 

with (^(m) a smooth function and we used a compact notation to denote the double sum 
over i and j. 

Again taking ip{m) = 1 we have that f{m,t) is a p.d.f. for all times, whereas the 
choice (/?(m) = m yields the evolution of the mean number of mutant cells. It is a simple 
computation to show that we get the same equation for the mean as in classical Luria- 
Delbriick setting since 

— I f{m,t)mdm = p{N{t),r]i) / p{m,'dj)['dj + r]i]f{m,t) dm 

ot Jr+ Jr+ 

= ^M{t)+nN{t). 

Let us now compute the evolution of the variance. Taking ip{m) = m? we have for the 
second moment M2{t) 



d_ 
di 



oo ,. 

f{m,t)m^dm = p{N{t),rji) p{m,{}j)f{m,t)['d]+T]f + 2rii^j 

+ 2m{r]i + dm, 

= <7f > +M{t){l3 + 2l3iiN{t) + 2^N{t)) 
+ /3(2 + /3) / f{m,t)m^dm. 



which gives the following expression for the kinetic variance 

= 2^V{t) + ^^M2{t) + /3M(t)(l - 2fiN{t)) + fiN{t){l + fiN{t)). (27) 

The above expression coincides with (jlip except for the presence of the additional term 
f3M{t) on the r.h.s. As before we omit for brevity the exact solution to (j27|) . 
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3.2 The generalized Fokker-Planck limit 

Following the analysis of section 12.21 we can introduce the mean field scaling defined by 
dH and (HH). 

Again the scaling does not affect the evolution of the time rescaled mean M(r) which 
in the asymptotic limit is governed again by ()15p . 

Now the time rescaled variance V{t) when e ^ is governed by 



dV{T) 
dT 



2jV{t) + -iM{t) + uNoe 



,71 



(28) 



The exact solution of the above equation for V{0) = yields the classical expression of 
the variance in the Lea-Coulson formulation ^1] 



V{t) 



71-7 



'^^(e(7l-27)T _ 1) 

71 - 27 



-'-fT 



7 

— A^oe^^fe'^^ - 1) - i^NoTe^\ 
I 7 



27 / 71, 7 7^ 71 
27 = 71, 

7 = 71- 



(29) 



Now let us consider the evolution of the rescaled probability density function /(m,r). 
We have 

d C 1 ^^'^ /* 

— / f(m,T)(p{m) dm = - y p(N(T),rji) / p(m,'d j)f(m,T)[ip(m') — ip(m)] dm. (30) 

By the same method introduced in section 12.21 we expand ip{m') in Taylor series and 
substitute into ([26]) to get 

— I f{m,T)(p{m)dm = -y^—'y'p{N{T),r]i) / p{m,^j)f{m,T){'dj+r]ifip^''\m)dm. 



k=l i,j=Q 

Using the binomial formula on (i?j + r]i)^ we can write 

00 ^ A; 



^ J fin^,rMn^)dm = ^l_^Q 
'''^ k=l h=o ^ ^ 



£ 



< 



>^ f{m,T)if^^\m)dm, 



where < > and < -rf^^^ > denote the /i-th moment and the {k — /i)-th moment of 
and r] respectively. 

We are interested in taking the limit e, /3i, /?, /x 0. Now since rn and "dj are given by 
Poisson processes 



< rj^-h >= ^n{t) + 0{n'^), < r >= (3m + 0{/3'^), 
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and then 

lim = < 7m, h = k; 

' [0, 0<h<k. 

In the hmit we finally obtain 

00 ^ ^ 

(m) dm 



d_ 



Jm+ ^ Jr+ 

°° 1 r 

"^i'^^yin / f{m,T)ip^''\m)dm, 

r-^ ^' Jr+ 



k=l 

which corresponds to the generalized Fokker-Planck equation 

d 



^^f{m,T) = £KM{{im + uN{T))f{m,T)), (31) 



where Ckm{') is the Kramers-Moyal operator given by ()2ip . 

Remark 1 (A simplified model) A further explicit analysis for equation \31\) seems 
quite difficult. Let us remark that taking the simplified microscopic dynamic 

m! = m + -d + T], (32) 

where now the growth of the random variable m depends only on the mean value of the 
variable itself, in weak-form we obtain the kinetic equation 

r- CO 

O 

di 



f{m,t)(p{m)dm = p{N{t),r]i)p{M{t),'dj) / f{m,t)[(f{m')-(f{m)]dm. (33) 



In the same scaling limit we obtain 

f f(rn,T)(f(rn)drn = — (i/N(t) +jM(t)] f f{m,T)Lp''^\m) dm, 



d_ 



which corresponds to the generalized Fokker-Planck equation 

— /(m,r) = (z.iV(T) +7M(r)) CKM{fim,T)). (34) 

Setting g{m, r) = f{m, r) and passing to the Fourier transform we obtain readily 

— 5(^,r) = (i/iV(r)+7M(r))5(e,r)(e-^«-l). 
The solution can be written in the form 

9{C, t) = m) exp (^(e-^« - 1) (^uN{z) + jM{z)) dz^ , (35) 
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which now can be integrated exactly and gives the characteristic function 

m, r) = m) exp ((e-^« - l)M(r)) . (36) 

Expression h30^ coincides with the characteristic function of a Poisson process with mean 
M(t). Note however that in our case m is a continuous variable and so we cannot conclude 
that the unknown mutant distribution has a Poisson distribution but it would be rather 
characterized by some continuous approximation of a Poisson distribution, like incomplete 
Gamma functions l2Tj. 

3.3 Diffusion approximations 

In order to compute the distribution of mutants we can clearly use the characteristic 
function (|24p or (j36p . However if we are interested in developing more realistic models 
where mutations are time dependent and related to space variables, like in a multicellular 
organism, we cannot rely on explicit solutions. Thus we must either solve the kinetic 
models d?]) and (|26p (or (j33p ). under the scaling (jl4p . or the corresponding generalized 
Fokker-Planck equations ()20p and ()3ip (or ()34|) ). In the latter case the Kramers- Moyal 
operator (j20p must be truncated at a finite number of terms and solved numerically. 

The Pawula theorem |27j states that truncations of the generalized Fokker-Planck 
equation in the form (j20p with finite derivatives greater than two leads to a contradiction 
to the positivity of the distribution function. It is then natural to consider a second order 
truncation to (|20p which in the original variables corresponds to the following diffusion 
approximation of the Luria-Delbriick dynamic 

§:^h^^^) + ^[il^ + ^^{T))f{m,T)] = ^uNir)^fim,r). (37) 

It is immediate to verify that the solution to the above equation has the same mean and 
variance as the original Luria-Delbriick distribution. 

Similarly we obtain the following diffusion approximation of the Lea-Coulson process 

^f{m,T) + -^[{jm + i^N{T))f{m,T)] = ^^[{jm + uN{T))f{m,T)]. (38) 

Clearly the above Fokker-Planck approximation preserve the mean and the variance evo- 
lution of the classical Lea-Coulson setting. 

Note that for the simplified model we obtain the standard diffusion approximation of 
a Poisson process 

—fim, t) + {^M{t) + vN{T)) — f{m, r) = -{^M{t) + z.iV(r))-^/(m, r). (39) 
OT am I om'^ 

We remark that higher order truncation can be considered as well. For a detailed 
discussion of the properties of solutions of Kramers-Moyal-expansions for a discrete Poisson 
process, truncated at an arbitrary order we refer to [28]. Although positivity of the solution 
is lost truncations at higher order are in better agreement with the exact solution than 
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the second order solution. However the numerical solution of higher order truncations 
presents some non trivial difficulties related to stability of the solution. For such reason 
here we will not explore further this direction. 

Let us finally remark that equations ([37|) . (j38|) and (p9|) are of the general type 

^/(m,t) + + = ^[(c(t)m + d(i))/(m,t)]. (40) 

In the case c{t) = they can be easily reduced to the heat equation introducing the 
change of variables (see for example [29] ) 



f{m,t) = u{t)g{y{t,m),T{t)), 



u{t) = exp ^— J a{s) ds 



y{t,m) = u{t)m — I b{s)u{s)ds, T{t) = / d{s)u{s)'^ ds , 
Jo Jo 



:9iy,r) = ^giy,T). (41) 



which gives 

d 

Thus from the well-known solution to the heat equation we recover the explicit solution 
for the probability density of mutants in time for equations (j37p and ()39p . We omit the 
details. 

From the above considerations it is clear that approximations (|37p . (j39p and (j38p may 
be inadequate since the distribution of mutants, in general, is far away from having a 
Gaussian shape. A possible alternative is to replace in the mutation dynamic the Pois- 
son process with a different stochastic process in such a way that the scaled mean-filed 
equations in the asymptotic limit originate Fokker-Planck models of the type studied in 
[71 [26]. Such Fokker-Planck equations cannot be reduced to the heat equation and origi- 
nate Gamma-like distributions in the long-time behavior. Here we do not explore further 
this direction leaving it open for future research. 



4 Numerical examples 

In this section we compare the continuous distribution of mutants obtained using the 
different kinetic models ([7]) and (I26p in the generalized Fokker-Planck limit and some 
standard methods to compute the approximated discrete distribution in the Luria-Delbrck 
and Lea-Coulson settings (see Lemma 2 page 11 in [31]). 

For the numerical solution of the kinetic master equations we adopt a standard Monte 
Carlo simulation method as it usually done in rarefied gas dynamics (see [6], 25j). Here 
the main difference is the necessity to generate Poisson samples at each time step. This 
can be easily achieved by standard algorithms, see for example [T71 [2T]. 

The test cases considered has been proposed in fBTj- We start from an initial conditions 
where no mutants are present fo{m) = (5(0) and A^o = 1- The parameters are ^ = 10"'', 
a — /Li = 3 and the final computation time is i = 6.7. To reduce fluctuations the total 
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number of simulation samples is 5 x 10^. First we consider the Luria-Delbriick case for 
/3 = 2.5 and then the Lea-Coulson case for f3 = 2.8. 

In Figure [21 we report the solution obtained simulating the kinetic model d?]) for 
different values of the scaling parameter e. The results show the convergence of the mutant 
distribution prescribed by model towards the classical Luria-Delbriick solution obtained 
using Lemma 2 in [31j . 

Similarly we report in Figure [3] the solution of the kinetic model (126p for different values 
of the scaling parameter e. As expected convergence of the mutant distribution towards 
the Lea-Coulson solution, computed using Lemma 2 together with numerical quadrature 
as in [31], is observed. 

5 Conclusions 

We have derived kinetic and mean field equations corresponding to two most famous math- 
ematical formulations of the Luria-Delbriick experiment. First we discussed the original 
Luria-Delbriick formulation and have shown that under a suitable mean field scaling the 
master equation yields a generalized Fokker-Planck equation having the Luria-Delbriick 
distribution as solution. Next we extended our approach to the Lea-Coulson case and 
derived the corresponding kinetic model and asymptotic generalized Fokker-Planck limit. 
Following the theoretical analysis we presented results of numerical experiments conducted 
using a Monte Carlo method. We have shown a computational evidence of convergence of 
our models towards the classical formulations. 

Beside the mathematical aspects of the problem, the interest in the continuous mod- 
eling of the mutation process presented here is twofold. On the one hand, the models 
provide a new environment for the development of numerical methods to compute mutant 
distribution, and, on the other hand, they represent a first step towards the development 
of more realistic mean field models where the distribution function depends also on the 
physical location of the mutant cells. 
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Figure 2: Luria-Delbriick setting. Distribution of mutant cells at t = 6.7 with /3 = 2.5 for 
e = 0.1 (top) and e = 0.01 (bottom) in the kinetic model ([7]). The reference solution is 
computed using Lemma 2 in [31] . 
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Figure 3: Lea-Coulson setting. Distribution of mutant cells at t = 6.7 with /? = 2.8 for 
e = 0.1 (top) and e = 0.01 (bottom) in the kinetic model (j26p . The reference solution is 
computed using Lemma 2 in [31] and numerical quadrature. 
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